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Abstract. It is well-known that the convergence of Krylov subspace methods to solve linear sys- 
tem depends on the spectrum of the coefficient matrix, moreover, it is widely accepted that for both 
symmetric and unsymmetric systems Krylov subspace methods will converge fast if the spectrum of 
the coefficient matrix is clustered. In this paper we investigate the spectrum of the system precon- 
ditioned by the deflation, coarse grid correction and adapted deflation precondit loners. Our analysis 
shows that the spectrum of the preconditioned system is highly impacted by the angle between the 
coarse grid space for the construction of the three preconditioners and the subspace spanned by the 
eigenvectors associated with the small eigenvalues of the coefficient matrix. Furthermore, we prove 
that with a certain restriction the accuracy of the inverse of projection matrix also impacts the 
spectrum of the preconditioned system. Numerical experiments emphasized the theoretical analysis. 

Key words, spectrum, coarse grid space, deflation, preconditioner, perturbation analysis, pro- 
jection matrix, iterative solvers, domain decomposition 



1. Introduction. We consider the iterative solution of the hnear system 

Ax = 6, 

where A E M"^" is symmetric positive definite(SPD). It is well known that the con- 
vergence of Krylov subspace method for solving linear system highly depends on the 
eigenvalue distribution. Recently many studies [l-3. il Onl 7j have shown that removing 
the very small eigenvalues makes the spectrum more clustered and consequently the 
convergence is improved. Preconditioners based on deflation technique and coarse 
grid correction are two main types of the spectral preconditioners that deflate or shift 
small eigenvalues. 

To describe the preconditioners discussed in the following sections, we define 
(1.1) E = Z^AZ, ZeM"^''. 



The matrix of the form ( 1.1 1 is referred to as projection matrix. The subspace spanned 
by the columns of Z is referred to as coarse grid space in domain decomposition 
method or deflation space in deflation method. In ideal case, the coarse grid space 
or deflation space should contain the vectors corresponding to the lower part of the 
spectrum that are responsible for the stagnation in convergence. Section 2 will show 
that the preconditioned systems have good properties when the preconditioners are 
constructed with the ideal subspace. For simplicity, we refer to the subspace spanned 
by the columns of Z as coarse grid space throughout the paper. In particular, we 
refer to the coarse grid space spanned by the eigenvectors associated with the small 
eigenvalues of A as the exact coarse grid space. 

In [21 [TS] , the deflation preconditioner is defined by 

(1.2) Pd = I - AZE-^Z^. 

Obviously PjjA is singular since Pd is singular. Fortunately, Krylov subspace methods 
can still solve singular linear system as long as it is consistent; furthermore, zero 
eigenvalues do not impact the convergence since the corresponding eigenvectors never 
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enter the iteration. To measure the convergence rate of Krylov subspace methods for 
symmetric positive semidefinite system, the effective condition number [TS], Keff, is 
defined by the ratio of the maximal to the minimal nonzero eigenvalues. 

In contrast to the deflation method, the preconditioners based on coarse grid 
correction shift the small eigenvalues to the large ones. In [3], a well-known coarse 
grid correction preconditioner in domain decomposition method is defined by 

(1.3) Pc = I + ZE-^Z'^. 

The abstract additive coarse-grid correction is M^^ + ZE^^Z'^ , where M is the sum 
of the local solves in each subdomain. Another well-known preconditioner in domain 
decomposition method is the abstract balancing preconditioner [14J 

(1.4) Pbnn = {I - ZE-^Z'^ A)M-'^{I - AZP-^Z'^) + ZE-^Z"^. 
The adapted deflation preconditioner [3j is defined by 

Padefi = M^^Pd + ZE-^Z'^. 

It has been shown in [3] that Padefi is less expensive than Pbnn but as robust as 
Pbnn- Moreover, Pbnn A and Padefi A have identical spectrum. In both Pbnn 
and Padefi, the first item is responsible for the fine grid and the second item term 
for the coarse grid. So they are called two-level preconditioners. In this paper, we 
confine our analysis on one-level methods. Thus we let M — I and define Pa by 

(1.5) Pa = I - AZE-^Z'^ + ZE-^Z^. 

Let X be a basis of coarse grid space. Obviously the following equation 

X{X^AX)-^X'^ ^ Z{Z^AZ)-^Z^ 

holds in exact arithmetic since there exists a nonsingular matrix C € W^"^ such 
that X — ZC. That means Pd, Pc and Pa are determined uniquely by coarse 
grid space. Thus we can choose an appropriate basis to form Z and then construct 
preconditioners. In section 3, we will show this fact is helpful for the analysis. 

The paper is organized as follows. In section 2, we briefly review the spectral 
properties of the preconditioned system when all preconditioners are constructed with 
the exact coarse grid space and inverse of projection matrix. Section 3 presents the 
perturbation analysis on the spectrum of the preconditioned system when the three 
preconditioners are constructed with the coarse grid space spanned by the approx- 
imate eigenvectors. Section 4 anlyses the spectrum of the preconditioned system 
when the inverse of the projection matrix are solved inexactly. Numerical results are 
presented in Section 5. 

2. Coarse grid subspace spanned by the exact eigenvectors. Let (Ai,Vi) 
be an eigenpair of A and Vi be normalized, (w^, • • • , «„) is orthogonal since A is SPD. 
The spectral decomposition of A can be written as 

(2.1) A^{V,V^) 

where A = diag{Xi, • • • , A^}, Aj^ = diag{Xr+i, • • • , A„}, V = (wi, • • • , Vr), and V± = 
{vr+i, ■ ■ • , Vn)- In this section, we assume that the eigenvalues of A are the small ones 
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to be removed and that the three preconditioners are constructed with Z — V and 
the exact inverse of projection matrix. 
Theorem 2.1. Let 

(2.2) Pb=I - AVE-^V'^ , E = V^AV. 

Then PdA is symmetric and has the following spectral decomposition 

PDA^iV,V^) 

In addition, PdA = APjj. 

Proof. Obviously, PdA is symmetric and PdAV = 0. From the definition of 
we have 

E = V^AV = V'^VK = A. 

Then we obtain 

PdAVi_ = AVi_ - AVK-^V'^AVa^ 

Thus the first result holds. Similarly, we have APpV = AV - AVKk-^V'^V = and 
APdV± = AVi_ + AVK-^^V^Vi_ = Vj_A_l. Thus PdA = APd- □ 
Theorem 2.2. Let 

(2.3) Pc ^ I + VE-'^V'^ , E = V'^AV. 

Then PqA is symmetric and has the following spectral decomposition 



In addition, PcA — APc ■ 

Proof. The proof is similar to that of Theorem 2.1 □ 
Theorem 2.3. Let 

(2.4) Pa^ I ~ AVE-^V'^ + VE-^V'^ , E^V^AV. 

Then PaA is symmetric and has the following spectral decomposition 

PaA^{V,Vi_) 

In addition, PaA = APa- 

Proof. The proof is similar to that of Theorem 2.1 □ 

Theorems 2.1|2.3 show that the small eigenvalues are removed no matter whether 
the preconditioners are applied on the left or right side of A, moreover, the rest 
eigenvalues and all eigenvectors are not changed. In addition, it follows from these 
theorems that Kef /{PdA) < k{PcA) and Kef f [PdA) < k{PaA), where k is the 
spectral condition number P]. 
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3. Coarse grid subspace spanned by the approximate eigenvectors. In 

practice for a large system, it is ratlier expensive to obtain tlie exact coarse grid 
space spanned by eigenvectors. An alternative to make the coarse grid space have 
rich information on eigenvectors is to build a large subspace[Tl O EH [H]- However, 
this way will lead to a large projection matrix whose inverse is also expensive to 
solve exactly. So it is natural to construct preconditioners by replacing the exact 
coarse grid space by an approximate one and computing E^^ approximately. We 
hope that the system preconditioned by the resulting preconditioners has the similar 
eigenvalue distribution as described in the previous section so that the convergence 
can be improved. That motivates us to analyse how the perturbation in the coarse 
grid space impacts the spectrum of the preconditioned system. 

Assume Z is column orthogonal in the following discussion. Let Z denote the 
subspace spanned by the columns of Z. Let Z-^ be the orthogonal complement of Z 
and Z± be an orthogonal basis of Z-^. Likewise, let V be the subspace spanned by 
the columns of V, be the orthogonal complement of V and V± be an orthogonal 
basis of . 

Let a denote the singular value of a matrix. Let dist{Z, V) denote the dis- 
tance between subspaces Z and V. It has been proven in [71 Theorem 2.6.1] that 
dist{Z,V) = amaxiZ'^V±) = (7ynax{V^ Z±). Note that < dist{Z,V) < 1 since Z 
and V have the same dimension. We define the acute angle between subspaces Z and 
V as 6* = aics'm dist{Z,V). 

Lemma 3.1. Let he the acute angle between subspaces Z and V that have the 
same dimension. Let Z and V be the orthogonal bases of Z and V respectively. Then 

sin 6* = (TmaxiZ'^V±) = a;nax{V'^ Z±), 
COS0 = a,mn{Z'^V) ^ (T,mn{ZlV±). 

Proof. The proof of the first equation can be found in [7, Theorem 2.6.1]. Since 
{V,V±) is orthogonal, it follows from ]| (y, Vj_)"'"Za;]|2 — 1 for all unit 2- norm a; e M'' 
that WV^ZxWl + \\VlZx\\l = 1. Thus 

(JrmniV^Zf^ miu Wv'^ Zx\\l = l~ max \\VlZx\\l 

|a;||2 = l ||2;j|2 = l 

= 1 - a„,ax{Vlzf ^cos^ 9. 

Similarly, since (Z, Zj^) is orthogonal, it follows from ]|(Z, Z±)'^V±x\\2 — 1 for all unit 
2-norm x G Rt""'') that WZ'^V^xWl + \\ZjV^x\\l = 1. Thus 

(t™„(ZIV1)2= min \\ZlV^x\\l = l- max \\Z'^V^x\\l 

||^||2— 1 ||a;||2 — 1 

= l~a,nax{Z'^V±^f =COs'e. 

Thus, COSO = armniZ^V) = (J,mn{ZlV±). □ 



Let Pd be defined by ( 1.2 1. PdA is an symmetric matrix since A is SPD. Hence, 
Courant-Fischer Minimax theorem[7] can be applied to estimate the eigenvalues of 
PdA, which is described in Theorem 13.21 



Theorem 3.2. If A e R"^" is symmetric, then 

\k{A) = max min x'^Ax, k — 1, 

dirn{S) — k S, || a; || 2 — 1 
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Theorem 3.3. Let Pjj be defined by (1.2) and 9 be the acute angle between V 
and Z. Then PjjA has r zero eigenvalues, moreover, if cos 9 ^ 0, then the nonzero 
eigenvalues of P^A satisfy 

Amm(A_L) - e_D < X(PdA) < \nax(-^±) + VD, 

whereijD = Amax(A_L)(sin e'+sin^ 6*) andeo = r]D + \\E^^\\2{\\E\\2+>^maxi^±)y ta.^'^ 9. 

Proof. It follows from PoAZ that PdA has r zero eigenvalues. From Theorem 
|2.H we have 

(3.1) x^PdAx = x'^V^KaYIx. 



For all unit 2- norm x e M", it can be written as x 
and X2 e . Moreover, there exist t € M"" and s G 

X2 = Zj_S. 



xi + X2, where xi & Z 
such that xi = Zt and 



x'^AZE-^Z'^Ax = xlAZE-^Z'^ Axi + x'^AZE-^Z'^Ax2 
+xlAZE^^Z'^Axi +xlAZE-'^Z'^ Ax2 
= t^Z^ AZt + t^Z'^ Ax2 + xlAZt 

+xlAZE-^Z'^Ax2 
= xfAxi + x^Ax2 + xT^Axi + xT^AZE^^Z'^ Ax2 
= x^Ax - X2AX2 + X2AZE^^Z^Ax2. 



Then we obtain 
(3.2) 



x'^PdAx = x'^ Ax 



x^AZE-^Z^Ax 
= xlAx2 - xlAZE-^Z'^ Ax2. 



Substract (3.11 from (3.2) on both sides, then we obtain 

(3.3) x'^PdAx = x'^PdAx + x'^Ax2 - x'^V^kiVlx ~ x'^ AZE^^ Z'^ Ax2. 



X2 Ax2 - Vj_A_lVj' X 



x^ Ax2 - {xi + xi )Vj_A_lVj' (.Ti + X2) 



- = X2AX2 — {xj + xj^ 
= x'^VAV'^X2 - FlAj^fJ Xi 

~xJVj_Ai_Vlx2 - x'^Vj_Ai_Vlxi 
= s'^{ZlV)A{V^Z^)s - t^{Z^V^)A^{VlZ)t 
~t'^iZ^V^)A^{VlZ^)s - s'^{ZlV^)A^{VlZ)t. 



W^ih = \\t\\2 and IIX2II2 
we obtain 



|s||2 since Z and Z± are orthogonal. Using Lemma 3.1 



(3.4) 



\x^Ax2 - x^V±A^Vlx\ < (|!x2||^||A||2 + ||a;i||2||Ai||2)sin2 9 

+ l|a;i||2||a;2||2||A_L||2sin26'. 



For the last item on the right side of (3.3), we only need to estimate the bound 
of ZlAZ because of x'^AZE-^Z^Ax2 = s^ZlAZ)E-\ZlAZ)'^s. It follows from 
(Z,^j_)(Z,Z_L)'^ = /that 



AZ = Z{Z'^AZ) + Z^{Z'iAZ). 
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Multiply Vf from left on both sides of above equation, then we have 

{VlZ^_){ZlAZ) = VlAZ - iVlZ){Z^AZ) 

= A±{VlZ) - {VlZ){Z^AZ). 

VfZ± is invertible since cos6' 7^ 0. Thus 

(3.5) (ZlAZ) = (VlZ^r'A^iVlZ) - (Vf Z^)-'{Vl Z)E. 



Using Lemma ISJ] we have \\ZlAZ\\2 < i\\E\\2 + ||A_l||2) tan0. Thus 

(3.6) x^AZE-^Z^Ax2 < ||a;2||^||-B-i||2(||^||2 + II A^lb)' tan^ 9. 
E^^ is SPD since A is SPD. Consequently, 

(3.7) xlAZE-^Z'^Ax2 = {Z'^ AX2Y E'^ {Z^ AX2) > 0- 



PdA is symmetric since A is SPD. Applying Theorem 3.2 to 3.3 with (3.4) and 



(3.7), we have 

\{PdA) < X{PdA) + \xlAx2 - x'^V^KiYlx] 

< ^niax{.J^±) + A„aa;(A_L)(sin6l + sin^ 9). 



Note that PdA has r zero eigenvalues. Applying Theorem 3.2 to 3.3 with (3.4 1 and 



(3.6), the lower bound of the nonzero eigenvalues of PdA is given by 

KPdA) > \{PdA) - \xlAx2 - x'^V±Kiylx\ ~ x^AZE-^Z'^Ax2 

max 

(Aj^)(sin6' + sin^ 9) 

'\\E-^U\\E\\2 + XmaA^±)fi''i^^0. 

Thus the theorem holds. □ 



Theorem 3.3 shows that in exact arithmetic, as 9 approaches to zero, the max- 
imal and minimal nonzero eigenvalues of PdA converge to XmaxiA±) and Xmin{A±) 
respectively. Hence with an appropriate coarse grid space the spectrum of PdA will 
be similar to that of PdA. Considering rounding error, however, PdAZ may be not 
equal to zero matrix, which implies PdA possibly has some eigenvalues around zero 
that should be equal to zero in exact arithmetic. So there is a risk for Pd to lead to 
a poor spectrum of the preconditioned system. We just mention rounding error here 
since roundoff analysis is outside the scope of this paper. 

The authors of [3 has shown that PdM^^A can be related to Padefi^ in terms 
of their spectrum, which is described in the following theorem. 

Theorem 3.4. Assume M is an SPD matrix. Let the spectrum of PdM~^A 
be given by {0, . . . , 0, 7r+i, ■ • • , 7n} with jr+i < 7r+2 < ■ • • < 7n- Let the spectrum 

of {M^^Pd + ZE^'^Z^)A fee {1, ... , 1,/ir+l, . . . with flr+l < Mr+2 < ■ • • < Mn- 

Then, 7^ — /i^ for all i = r + 1, . . . ,n. 

Proof. The proof can be found in [31 Theorem 3.3]. □ 

Corollary 3.5. Let the .spectrum of PdA be given by {0, . . . , 0, A^+i, . . . , A„} 
with Xr+i < Xr+2 ^ • • • ^ Xn- Then the .spectrum of Pa A is {1, . . . , 1, A^+i, . . . , A„}. 
Proof Let M = I. The result follows from Theorem □ 



Theorem 3.6. Let Pa be defined by (1.5). Let 9 the acute angle between subspaces 
Z and V. If con 9 7^ 0, then the eigenvalues of Pa A satisfy 

min{l, Atoj„(A_l) - ed} < X{PaA) < max{l, Amai;(A_L) +ryD}, 
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where rjo and en are defined in Theorem p7 

Proof. It follows from the combination of Theorem |3.3| and Corollary |3.5| □ 
PcA is not necessarily symmetric although both A and Pc are symmetric. So we 

can not apply Courant-Fischer Minimax theorem to estimate the eigenvalue of PcA. 

Fortunately, the eigenvalues of PcA are positive since Pq'A is similar to a SPD matrix. 

Thus the eigenvalues of PcA can be bounded by max{x'^PcAx} and min {x'^PcAx} 

for all unit 2-norm a; G K". 

Theorem 3.7. Let Pc be defined by (1.3). Let 6 be the acute angle between V 

andZ. If cosd^O, then 



KnaxiPcA) < max{l + A,; 
^min{Pc^) > min{l + A 



.X 

(A), A 

max (A_l)} + ec, 
(A),A™„(A_l)} - ec, 



where ec = \{\raax{J^i)\\E^^\\2 + 1) tan 6* + sin 6* + sin^ 6'. 

Proof PcA is similar to A + A^^'^ZE^^Z^ A^l"^ since A is SPD. Moreover A + 
A^I'^ZE-^Z'^ A^l'^ is SPD as well since A and E arc SPD. Thus the eigenvalues of 
PcA are positive. 

For all unit 2-norm x G M", it can be written as a; = xi + X2, where x^ ^ Z 
and X2 G Z^. Moreover, there exist t e M'" and s e M"^'" such that x\ = Zt and 
X2 = Zj_s. 



(3.8) 



x'^PcAx = x'^Ax + {xi + X2f ZE-^Z'^ A{xi + 2:2) 



'^Ax+{xi +X2fz 

= x'^Ax + x'^ZE-'^Z'^Axi + x^ZE-^Z'^ Ax2 
= x^ Ax + x^xi + x^ZE^^Z^ Ax2. 



From the definition of Pc in (|2.3|), we obtain 
(3.9) 



iFPcAx = x'^Ax 



x^VE-W^Ax 
= x^yla; + x'^VV'^x. 



Subtract (|3.9|) from (|3.8|) on both sides, then we have 
(3.10) 



x'^PcAx 



x'^PcAx 



xlZE^'^Z'^ Ax2 + xi - x^yi/^a; 



= xlZE-'^Z'^ Ax2 + KL^f xi 
-a;f yy'^X2 - 2xlVV'^X2. 

Since both A and E^^ are symmetric. 



Z£:-iz'^Aa;2 = x^AZE-^Z^^x^ 



s^ZlAZE- 



Using (3.5), we have 

{ZlAZ)E-^ = {VlZ^)-^k^{VlZ)E-^ {VlZ^rWlZ). 



Note that ||a;i||2 = \\t\\2 and ||a;2||2 = ||s||2- Using Lemma 3.1 we obtain 

(3.11) \xlZE-^Z'^Ax2\ < ||xi||2||a;2||2(||Ai||2||i;"'||2 + l)tan0. 
Since ||a;||2 = 1, we have the following bounds with Lemma 3.1 

0<x'^V±Vlxi < ||xi||^sin2 6', 

(3.12) 0<x^VV^X2<\\x2\\lsm^e, 

\xJVV^X2\ < ||a;i||2||x2||2sin0. 
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With ( [3I0I ), ^Jl\ and ^J^, we obtain 



XrmniPcA) ~ Cc < Pc Ax < XrnaxiPcA) + Ec, 

where ec — ^{Xmax{-^±)\\E^^\\2 + 1) tanfl + sin 6* + sin^ 6. Thus the theorem follows 
from min{a:'^Pc^a;} < A(Pc^) < T[iayi{x'^ Pc Ax} for all unit 2-norm x S E". □ 



Theorem |3.7| shows that the maximal and minimal eigenvalues of PcA converge 
to max{Amaa;(A_L), 1 + XmaxW} and min{A™„(A_L), 1 + Amin(A)} respectively as 9 
approaches to zero. Therefore it can be expected that with an appropriate coarse grid 
space PcA will have the similar spectrum of Pc'A. Analogously, it can be proved that 
the maximal and minimal eigenvalues of APc have the same bounds as that of PcA. 

4. Inexact inverse of projection matrix. For a large projection matrix, it is 
expensive to solve its inverse exactly. So we want to replace a projection matrix by 
an approximate one that is cheap to compute its invers e. I n this seciton, we assume 
H is an invertible matrix as an approximation to E in (2.2 1. 

Theorem 4.1. Let Pd = I - AVH^^V'^ and p = EH^^ - I. If the eigenvalues 
of PjjA are real, then 

-^D < HPdA) < A 

max 

w/iere ^r, = IIPII2IIAII2. 

Proof. \\p\\2 measures the distance between H and E when H is used as the right 
inverse of E. Obviousl y ||p ||2 = if and only \i E = H . Assume ||a;|j2 = 1 and use 
the definition of in (2.2 1, then we obtain 



jFPdAx = x'^Ax 
= x'^Ax 



x'^AVH^W^Ax 
x'^VV'^Ax - x^VpV'^Ax 



= x^PdAx - x^VpAV^x. 
Since \x^VpAV^x\ < WphWMh, 

-II/0II2IIAII2 < x^PdAx < A 

max 

{PDA) + \\ph\\Ah. 

Therefore the theorem follows from min {x'^PdAx} < X{PdA) < max{x^PD^a;}. □ 

Theorem 4.1 shows that if the eigenvalues of PdA are real, the maximal and min- 
imal eigenvalues of PdA converge to XmaxiA±) and respectively as ||p||2 approaches 
zero. That implies that PdA will have small eigenvalues around zero. As a result, 
the spectrum of PdA will become worse than that of PdA. 

_Theorem 4.2. Let Pc = I + VH'^V'^ and p ^ H^^E - I. If the eigenvalues 
of PcA are real, then 



Xmax{PcA) < max{l + A 

n 

XminiPcA) > min{l + A„ 



■ X 

(A), A 

max 

,(A),A™„(A_l)}-Cc, 



where ^C ^ \\ph- 

Proof. \\p\\2 measures the distance between H and E when H is used as the left 
inverse of E. Obviousl y \\p \\?. = if and only ii E = H. Assume ||a;|j2 = 1 and use 
the definition of Pc in (2.3 1, then we have 



c'^PcAx = x^ Ax 



x'^VH^^V^Ax 
x^Ax + x'^VE'^V'^Ax + x'^VpV'^x 
x^PcAx + x'^VpV'^x. 
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Since \x^VpV'^x\ < \\p\\2, 

{PcA) - ||p||2 < x'^PcAx < X 

max {PcA) + \\ph. 

Therefore the theorem foUows from min {x'^PcAx} < XiPcA) < ma,x{x'^PcAx}. □ 



Theorem 4.2 shows that if the eigenvalues of PcA are real, as p approaches to 
zero, the maximal and minimal eigenvalues of Pc'A converge to max{Amax(Aj^), 1 + 
Amaa:(A)} and miujAmin (A_l) , l + Amm(A)} respectively. That implies PcA and PcA 
have the similar eigenvalue distribution if ||p||2 is small enough. 

Theorem 4.3. Let Pa = I - AVH^^V^ + VR-^V'^. Let pi = ER-^ - I and 
P2 = H^^E — I. If the eigenvalues of Pa A are real, then 

Xmax{PAA) < max{l,A 

max 

X„iin{PAA) > min{l, Ami„(A_L)} - ^a, 

w/iere = II Pi 11 2 II All 2 + II P2 II 2. 

Proof. Assume ||x||2 = 1 and use the definition of Pa in (2.4), then we have 

x^PaAx = x^Ax - x'^AVH-W^Ax + x'^VH-W^Ax 
= x^Ax - x'^VpiAV'^x + x^Vp2V'^x 
-x'^AVE^W^Ax + x'^VE-W^Ax 
= x^PaAx - x^VpiAV^x + x'^Vp2V'^x. 

Since \x^VpiAV^x ~ x'^Vp2V^x\ < ||pi||2||A||2 + IIP2II2, 

X„,,n{PAA) - IIP1II2IIAII2 - IIP2II2 < x'^PaAx < X,r^ax{PAA) + ||pi||2||A||2 + ||p2||2- 

Thus the theorem follows from min{x'^ PaAx} < X{PaA) < max{x^ Pa Ax} . □ 



Similarly, Theorem 4.3 shows that the maximal and minimal eigenvalues of PaA 
converge to max{l, Amaa:(Aj_)} and min{l, Amin(A_L)} respectively as both ||pi||2 and 
IIP2II2 approach to zero, if the eigenvalues of PaA are real. That imphes PcA and 
PcA have the similar eigenvalue distribution if if is a good approximation to E. 

5. Numerical experiments. In this section, numerical comparison of various 
preconditioners will be reported. All tests were performed by using Matlab(R2010b) 
on an Intel Core2 Duo E7500, 2.93GIIz processor with 4Gb memory. Except for 
the deflation preconditioner, the system preconditioned by the coarse grid correction 
and adapted deflation preconditioners are not necessarily symmetric although A is 
SPD. Therefore we apply GMRES^ to solve the preconditioned system iteratively. 
In addition, we apply Gram-Schmidt method with reorthogonalization to maintain 
the orthogonality of basis of Krylov subspace[7| H]- 

5.1. Diagonal matrix. The first test case is a diagonal matrix with entries 10"'^, 
10"^ • • •, 10"\ 1, 10, 10.1, 10.2, • • •, 209.1. The right hand side is a vector of all 
ones. The matrix has 7 small eigenvalues less than 1 to be removed. The eigenvectors 
associated with these eigenvalues are the unit vectors, i.e., V = {ei,e2, ■ ■ ■ , 67), where 
Ci is the ith column of identity matrix. All tests are required to reduce the relative 
residual below 10^^^. The initial guess vector is always zero vector. GMRES method 
converges within 273 iterations without preconditioner. The perturbations in coarse 
grid space and projection matrix are generated by Matlab function rand. 

Table [O] shows the distance between the exact coarse grid space and the coarse 
grid space with various perturbation. It should be noted that it is expensive and 
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Table 5.1 

The distance between span{V} and span{V + rand/e} . 





e = le + Ol 


e = le + 02 


e = le + 03 


e = le + 04 


e = le + 05 


sin 6 


9.77e-01 


5.06e-01 


6.03e-02 


6.05e-03 


6.02e-04 




7.08e+01 


5.54e+01 


7.33 


5.78e-01 


4.43e-02 



unnecessary in practice to compute sin0 by Lemma 3.1 for a general linear system. 
Assume that (A^, Vi){i — 1, • • • , r) are Ritz pairs that are extracted from the perturbed 
coarse grid space by Rayleigh-Ritz procedure[8]. In general, maxi=i^...^r{\\Avi—XiVi\\2} 
denoted by reSmax in Table [5?T] decreases as the two subspaces approach to each other. 
So we can use reSmax to measure the distance between the two subspaces since it is 
cheap to obtain. 



Table 5.2 

The number of iterations of GMRES with various preconditioners that are constructed with 
Z = V + rand/e and E^^ . 





V, E-^ 


e= le + Ol 


e = le + 02 


e = le + 03 


e = le + 04 


e = le + 05 


Pd 


71 


>300 


>300 


>300 


109 


88 


Pc 


104 


273 


267 


222 


173 


144 


Pa 


72 


290 


231 


167 


117 


96 



10 10 10 10 10' 10 

eig(P A) 



10 10 10 10 10 10 

eig(P A) 




Fig. 5.1. The eigenvalue distribution of FiG. 5.2. The eigenvalue distribution of 

the system with the three preconditioners con- the system with the three preconditioners con- 
structed with Z = V + rand/ le + 03 and . structed with Z = V -i- rand/ le -{- 05 and E~^ . 



In Table |5.2[ the second column shows the number of iterations needed for con- 
vergence for GMRES with the three preconditioners in the case that both coarse grid 
space and i?^^ are exact. The columns 3-7 show the number of iterations required 
for convergence in the case that only coarse grid space is perturbed. As is shown, all 
preconditioners suffer from the perturbation in the coarse grid space if it is large(see 
the third column). As the perturbation decreases, Pc and Pa behave better, whereas 
Pd behaves better only if the perturbation is small enough. On the other hand, if the 
perturbed coarse grid space is close enough to the exact one(see the last two columns 



in Table 5.2 1, is more efficient than Pa and both of them are much more efficient 
than Pd- 
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Figure |5.1| and Figure |5.2| show the eigenvalue distribution of the system pre- 
conditioned by the three preconditioners. As discussed in previous section, because 
of rounding error, PdA has some tiny eigenvahies around zero. For this test case, 
rounding error does not impact the spectrum of PdA when the perturbation is smah 
enough. 

Table 5.3 

The number of iterations of GMRES with different preconditioners, where only projection ma- 
trix is perturbed and the perturbation of E is H = E + rand/e. 





e = le + 10 


e = le + 12 


e = le + 14 


e = le + 16 


Pd 


>300 


>300 


>300 


>300 


Pc 


111 


104 


104 


104 


Pa 


88 


88 


80 


80 


\\H-^E~I\\2 


1.68e-03 


1.08e-05 


1.77e-07 


1.23e-09 


\\EH-^-I\\2 


1.68e-03 


1.08e-05 


1.77e-07 


1.23e-09 



In Table |5.3[ the first three rows show the number of iterations required for con- 
vergence in the case that only projection matrix is perturbed. The distance between 
and are reported on the last two rows. As is shown, Pc and Pa are robust 
on the perturbation in projection matrix for all four cases, while Pu is sensitive to 
that even though is very close to E~^{see the last column). 



10 10 

eig(P A) 



Fig. 5.3. The eigenvalue distribution of FiG. 5.4. The eigenvalue distribution of 

the system with different preconditioners con- the system with different preconditioners con- 
structed with V and H = E -\- rand/le -\- 12. structed with V and = E + rand/le + 16. 



Figure 5.3 and Figure 5.4 show Pc and Pa successfully shift the small eigenvalues 
of A to the large ones; however, Po fails to deflate the small eigenvalues and generates 
much smaller ones, which leads to a worse eigenvalue distribution. 

5.2. Boundary value problem. We solve the following model problem 

-V • (kVw) = / in n = [0,l]^ 

u = on dfl, 



by the two-level multiplicative Schwarz method^2 . The model problem is discretized 
by FreeFem-|--|-[TBj and we obtain a coefficient matrix with the order of 10201. Tests 
are performed on irregular overlapping decompositions with the overlap of 2 elements. 
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These overlap decompositions are built by adding the immediate neighboring vertices 
to non-overlapping subdomain obtained by Metis |T5]. 

In the two-level multiplicative Schwarz method, one-level preconditioner is the 
restricted additive Schwarz preconditioner(RAS) 5j that is responsible to remove high 
frequency models of the original system, and the deflation, coarse grid correction and 
adapted deflation preconditioners are applied as two-level preconditioners that are 
responsible to remove lower frequency ones of the system preconditioned by one-level 
preconditioner. 

We choose Ritz vectors to span the coarse grid space, which are extracted from 
Krylov subspace during the solve of the system preconditioned by RAS. These vectors 
are approximations to the eigenvectors corresponding to the lower part of the spectrum 
of the system preconditioned by RAS. To enrich the information on lower part of the 
spectrum, we use Ritz vectors in a splitting way to construct the coarse grid space, 
see [IH3] and references therein. More precisely, let 



V 



Vll 

V21 



V12 
V22 



V2,r 
Vnparts,r 



store Ritz vectors columnwise, where nparts is the number of subdomains and r the 
number of Ritz vectors; let Zi store the orthogonal vectors obtained by orthogonalizing 
(vii, Vi2, ■ ■ ■ , Wi,r)- Then Z is formed in the following way 



Z : 









Z2 







''nparts 



Obviously, span{T^} is a subspace of span{Z}, moreover span{Z} is nparts times 
as large as spanjV}. So we can expect spanjZ} has richer information corresponding 
to small eigenvalues. Note that Z has a very sparse structure, which leads to a sparse 



structure of £^ in (1.1) as well 



Comparison of the three preconditioners are performed on two different configura- 
tions with highly heterogeneous viscosity. See [1] for details. Two cases are described 
as following: 

• skyscraper viscosity: for x and y such that for [9x]=0(mod 2) and [9y]=0(mod 
2), K = 10^([9?;] -I- 1); and k = 1 elsewhere. See Figure 5.5 

• continuous viscosity: K{x,y) = 10^/3 sin(47r(x -t- y) -I- 0.1). See Figure 5.6 



Figures 5.7||5.10 plot the convergence curves for the case of skyscraper viscosity 
against the various number of subdomains. As is shown, Pd has a long plateau in 
the convergence on the decomposition with 16 subdomains, while Pa and Pc improve 
convergence significantly. On the other decompositions, Pjj has almost the same 
number of iterations as Pa although the initial residual of Pjj is considerably less 
than Pa, moreover P^ and Pa performs better than Pc- In addition, two- level 
method with RAS and Pa and Pc outperforms one-level method with RAS on all 
four decompositions. 

Figures |5.11|5.14] plot the convergence curves for the case of continuous viscosity 
against the various number of subdomains. Again two-level method with RAS and 
the three preconditioners all outperform one- level method with only RAS. For this 
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Fig. 5.5. Skyscraper case 



10^ 




20 40 60 80 100 120 140 160 180 200 



Fig. 5.7. Skyscraper case: 16 subdomains, 
15 Ritz vectors spanning coarse grid space. 





— ■-. RAS 

1 RAS+P^ 
RAS+P(, . 




^ RAS+Pp 


\ \ 




\ \ 





20 40 60 80 100 120 140 160 ISO 



Fig. 5.9. Skyscraper case: 64 subdomains, 
20 Ritz vectors spanning coarse grid space. 
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Fig. 5.6. Continuous case 
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1 RAS+P^ 




__ RAS+Pj. 




^ RAS+Pjj 






\ \ 
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20 40 60 80 100 120 140 160 



Fig. 5.8. Skyscraper case: 32 subdomains, 
16 Ritz vectors spanning coarse grid space. 




20 40 60 80 100 120 140 160 180 200 



Fig. 5.10. Skyscraper case: 128 subdo- 
mains, 20 Ritz vectors spanning coarse grid 
space. 



case, Pd and P4 behave almost the same on all four decompositions and outperform 
Pc- 

Table |5.4| shows the maximal residual of Ritz pairs for both cases that are ex- 
tracted from Krylov subspace when solving the system preconditioned by RAS. 
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Fig. 5.11. Continuous case: 16 subdo- 
mains, 16 Ritz vectors spanning coarse grid 
space. 




20 40 60 aO 100 120 



Fig. 5.12. Continuous case: 32 subdo- 
mains, 20 Ritz vectors spanning coarse grid 
space. 




20 40 60 aO 100 120 140 



Fig. 5.13. Continuous case: 64 subdo- FiG. 5.14. Continuous case: 128 subdo- 

mains, 20 Ritz vectors spanning coarse grid mainss, 20 Ritz vectors spanning coarse grid 
space. space. 



Table 5.4 

Maximal residual of Ritz pairs for skyscraper and continuous cases. 



Nparts 


16 


32 


64 


128 


skyscraper 


1.152184e-06 


1.973786e-04 


2.417277e-02 


5.017373e-02 


continuous 


4.760946e-03 


2.395625e-02 


2.003093e-02 


2.494041e-02 



Note that the size of E becomes large when decomposition has 64 or 128 sub- 
domains. As a consequence, computing is costly and undermines gains in the 
number of iterations. Hence, we attempt to compute the inverse of an approximation 
to E instead of E^^ . Assume L and U are factors of incomplete LU factorization(ILU) 
with no fill-in of E. Note that E has a sparse structure because of the sparse structure 
of Z. L and U are also sparse and cheap to be inverted. Thus we replace E~^ by 
{LU)-\ 

Figure |5.15| and Figure |5.16| plot the convergence curves on decompositions with 
64 and 128 subdomains for skyscraper case. Pa and Pc are insensitive to the per- 



turbation in E ^, while Pd fails to converge. For the continuous case, Figure 5.17 



and Figure 5.18 show Pa and Pc behave stable on decompositions with 64 and 128 
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subdomains, while Pjj is only stable on decomposition with 64 subdomains since LU 
is almost equal to E. Table [53] and Table [5^ present the distance between LU and 
E for both cases respectively. 




10 20 30 40 50 60 70 30 90 100 50 100 150 300 250 300 



Fig. .5.15. Comparison of different projec- 
tion matrix on the decomposition with 64 subdo- 
mains for skyscraper case. 



Fig. 5.16. Comparison of different projec- 
tion matrix on the decomposition with 128 sub- 
domains for skyscraper case. 



Table 5.5 

The distance between LU and E for skyscraper case. 



Nparts 


64 


128 


\\E{LU)-' Ih 


3.8588e-08 


8.9712e+02 


\\{LUy'E-I\\2 


4.1433e-10 


8.6292 




Fig. 5.17. Comparison of different projec- FiG. 5.18. Comparison of different projec- 

tion matrix on the decomposition with 64 subdo- tion matrix on the decomposition with 128 sub- 
mains for continuous case. domains for continuous case. 



6. Conclusion. In this paper we present a perturbation analysis on the defla- 
tion, coarse grid correction and adapted deflation preconditioners when the inexact 
coarse grid space and inverse of projection matrix are applied for the construction of 
the preconditioners. Our analysis shows that in exact arithmetic the spectrum of the 
system preconditioned by the three preconditioners is impacted by the angle between 
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Table 5.6 

The distance between LU and E for continuous case. 



Nparts 


64 


128 


\\E{LU)-'-Ih 


4.1008e-14 


5.4090e-01 


\\{LU)-'E-Ih 


1.3890e-15 


1.0215e-01 



the exact coarse grid space and the perturbed one. Moreover, wc prove that with a 
certain restriction the coarse grid correction and the adapted deflation preconditioners 
are insensitive to the perturbation in projection matrix, whereas the deflation pre- 
conditioner is sensitive. Nmnc;ric;al results of the different test cases emphasized the 
perturbation analysis. For unsymmetric linear system, similar perturbation analysis 
on the coarse grid correction and adapted deflation preconditioners can be developed 
and requires further investigation. 
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